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We show that finite-range ahernatives to the standard long-range BKS pair potential for silica 
might be used in molecular dynamics simulations. We study two such models that can be efliciently 
simulated since no Ewald summation is required. We first consider the Wolf method, where the 
Coulomb interactions are truncated at a cutoff distance rc such that the requirement of charge 
neutrality holds. Various static and dynamic quantities are computed and compared to results 
from simulations using Ewald summations. We find very good agreement for rc ~ 10 A. For lower 
values of r^, the long-range structure is affected which is accompanied by a slight acceleration of 
dynamic properties. In a second approach, the Coulomb interaction is replaced by an effective 
Yukawa interaction with two new parameters determined by a force fitting procedure. The same 
trend as for the Wolf method is seen. However, slightly larger cutoffs have to be used in order to 
obtain the same accuracy with respect to static and dynamic quantities as for the Wolf method. 

PACS numbers: 02.70.Ns, 61.20.Lc, 61.20.Ja, 64.70. Pf 



I. INTRODUCTION 

Silica (8102) is the prototype of a glassformcr that 
exhibits a tetrahedral network structure. It is the ba- 
sic oxidic component of many minerals and technolog- 
ical glasses. In recent years, molecular dynamics (MD) 
computer simulations have provided valuable insight into 
static and dynamic properties of amorphous and crys- 
talline silicaJJ, Jl M, 11 [il [H, [H, Q 111, 
m, [O, E [ii Sllirill,!!! (mU^ in a classical MD 
simulation, the interactions between the atoms are de- 
scribed by an effective potential where, different from ah 
initio approaches, the electronic degrees of freedom are 
not taken into account explicitly. In many of the afore- 
mentioned MD studies, the so-called BKS potential [2^ 
has been used. Although it is a standard pair potential, 
it yields good agreement with experimental data. How- 
ever, it contains a long-range Coulomb interaction term 
and thus the computation of energies and forces is very 
expensive, in particular for large systems. Therefore, it 
is important to know whether silica can also be modeled 
by a potential with a finite range. 

The classical approach to evaluate Coulomb ener- 
gies and forces in a simulation is the Ewald summa- 
tion method where the sum over all Coulomb interac- 
tions is decomposed into a real space and a Fourier part 
[27l . m, [2§| . For a three-dimensional N particle system 
with periodic boundary conditions in all three spatial 
directions, Ewald method yields a N^/"^ scaling for the 
computational load [1^, [30| . Even worse is the case of 
a quasi-two-dimensional slab geometry for which Ewald 
method exhibits a iV^ scaling [sj. There are modifi- 
cations of the Ewald method such as particle-particle 
particle-mesh methods (PPPM) [s^,!!!] that yield a scal- 



ing proportional to iVlogA^, but these methods cannot be 
implemented as efficiently as in the three-dimensional 
case for the quasi-two-dimensional geometry [s^l- A 
scaling better than 0{N\ogN) can be reached by means 
of multipole methods [s^. However, these methods in- 
troduce a large computational overhead such that their 
use is only reasonable for a very large number of particles, 
say N > 10^ 

Apart from the case of quasi-two-dimensional geome- 
tries, there are other applications where it is difficult to 
handle long-range Coulomb interactions. For instance, 
the calculation of transport coefficients such as the shear 
viscosity via Green-Kubo relations is in general more 
complicated and less efficient when Ewald sums have to 
be considered [2^. This is also the case for Monte Carlo 
(MC) simulations. When using Ewald summation (or 
also PPPM) in MC simulations, the full Fourier part of 
the Coulomb energy must be evaluated after each parti- 
cle displacement by summing over N terms and Nk dif- 
ferent wavevectors used in the Ewald summation. Note 
that this is much worse than for MD where the sum over 
wavevectors can at least be used to compute the force 
over all N particles, while the energy must be recom- 
puted after each single move in MC simulations. 

For all these reasons, a reliable finite range pair poten- 
tial to simulate silica is highly desirable. In this paper, we 
address this issue by making the assumption that, due to 
screening effects, the long-range Coulomb interactions in 
typical ionic systems can in fact be truncated [37| . How- 
ever, a cutoff distance rc for a r^^ potential cannot be 
introduced in a similar manner as for a short-range po- 
tential since a truncated spherical summation over pair- 
wise Coulomb terms leads to a violation of charge 
neutrality. In order to circumvent this problem. Wolf et 
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al. [371 . |38[ have introduced a simple correction term in 
the truncated Coulomb sums that recovers the require- 
ment of charge neutrality. The Wolf method [13, S, [3§| 
has been used in recent simulation studies, e.g. for quasi- 
two-dimensional geometries [i^ , in MC simulations [4l[ , 
or for dipolar fluids [43 |. Of course, it is a priori not 
clear how large the cutoff radius Tc has to be in order 
to correctly reproduce the static and dynamic properties 
of the original model with long-range Coulomb interac- 
tions. We have carefully studied this issue using MD sim- 
ulations of amorphous silica based on the BKS potential. 
We show that a cutoff of about 10 A is necessary to ob- 
tain good agreement with the initial long-range model on 
a quantitative level. This indicates that, due to screen- 
ing effects, amorphous silica can indeed be described by 
an effective potential of finite range. Therefore, in our 
second approach we reparamctrizc the Coulomb term of 
the BKS potential by replacing it by a screened Coulomb 
(Yukawa) potential with two new parameters. The pa- 
rameterization of this Yukawa potential is achieved by a 
force fitting procedure based on previous MD simulations 
of BKS silica [3| . This alternative procedure is of course 
more involved than the simpler Wolf truncation, as the 
new pair potential must be carefully parameterized. 

The paper is organized as follows. In Sec. |TT] we in- 
troduce the truncated and screened Coulomb potentials 
studied in this work, and we describe how the free pa- 
rameters are fixed. In Sec. IIIII we use a given set of fixed 
parameters and compare in detail static and dynamic 
properties of the original long-range BKS model, and 
the two finite-range alternatives suggested in this paper. 
We present our conclusions in Sec. IIVI 



II. TWO FINITE-RANGE ALTERNATIVES 



In numerical simulations using periodic boundary con- 
ditions, the Coulombic part of the potential energy, i^coui, 
is given by the following formula: 
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Here, r^j 



Ti — Vj is the distance vector between particle 
i and particle j. The sum over n takes into account the 
interactions of a particle i with all the replicated image 
particles (including the images of particle i). Moreover, it 
is only conditionally convergent, i.e. the value of the sum 
depends on the order by which the terms are summed 
up. An efficient method to circumvent these problems is 
provided by the Ewald summation technique. However, 
finite range alternatives to Ewald sums are very desirable, 
as outlined in the Introduction. The analysis of such 
alternatives is the principal aim of this paper. 

We shall analyze the ability of finite-range potentials 
to reproduce the behavior of the original BKS model 
by comparing the results of molecular dynamics simu- 
lations of the new potentials against results employing 
Ewald sums. The latter results are taken from previ- 
ously published work [1, and they will be labelled 
"Ewald simulations" hereafter. Most of the Ewald sim- 
ulations were performed with N — 1008 particles, at the 
density 2.37g/cmf^ with Ewald parameters [a and N}~) 
optimized as in In particular the real space part 
of the Coulomb force is truncated at 10.17 A, and we 
take the results of these simulations as representative of 
the behavior of the "real" Coulomb potential. We know, 
however, that finite size effects are present for this system 
size, and we borrow additional data from the TV = 8016 
particles simulations from Q when necessary. 



A. The original long-range BKS model 

The functional form of the BKS potential is [1^1 

'/'of W =9a<Z0eVc(r)+A„^exp(-i?„0r)-^, (1) 

where a, (3 G [81,0], r is the distance between the 
ions of type a and (3. The values of the constants 
qa.qp, Aap, Ba/3, and Cap can be found in Ref. [1^. For 
the sake of computational efficiency the short range part 
of the potential was truncated and shifted at 5.5 A 
This truncation also has the benefit of improving the 
agreement between simulation and experiment with re- 
spect to the density of the amorphous glass at low tem- 
peratures [1, [!]■ Finally, the function 

Vc{r) = - (2) 
r 

is the Coulomb long-range term which will be further 
approximated by finite range potentials. The other terms 
in Eq. ([T|) will be unchanged. 



B. Truncation using the Wolf method 

As proposed in Refs. [33,[H,[3§|, the Vc{r) term in the 
BKS potential ([T]) may be approximated by the following 
form: 




while Vwir > rc) = 0. The potential ([U is a finite-range 
potential: Only particles separated by distances smaller 
than ?'c interact. Such a potential becomes computation- 
ally extremely useful when system sizes that are much 
larger than Tc arc studied. But Wolf's main point is that 
reasonable values of rc can lead to numerically accurate 
results (btI . That this is by no means a trivial statement 
can be appreciated in Fig. [T] which compares the initial 
Coulomb interaction in (jT]) to the expression given by 
Eq. ([4]). Both potentials are very different, and the trun- 
cation should in principle have a drastic effect. Indeed 
previous naive truncation attempts have been shown to 
produce quantitatively inaccurate results [2^. We shall 
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FIG. 1: Comparison of the pure 1/r Coulomb interaction, 
the Wolf truncation using Eq. Q for Tc — 10.17 A, and the 
Yukawa potential in Eq. (0, for Dy = 1.07, A = 5.649 A, 
rc = 10.17 A. 



demonstrate, however, that the Wolf method performs 
well in the case of liquid silica. 

The truncated form of the potential ^ can be justified 
as follows. Wolf [13 showed that the main error made 
when imposing naively a finite distance cutoff to a 1/r 
potential stems from the fact that a sphere of radius Vc 
centered around a particle is in general charged, so that 
its interaction with the rest of the system is not negli- 
gible. He proposed to approximate this interaction by 
using the potential 



Vir) 



1 



r < rc, 



which amounts to screening completely the charge con- 
tained in the sphere by placing the opposite charge on 
its surface. This form is also natural from the computa- 
tional point of view: Truncated potentials in molecular 
dynamics simulations are always shifted to avoid energy 
discontinuities at the cutoff, V{r = rc) = 0. Although 
the shift is sufficient to compute the energy, which was 
the initial problem considered by Wolf [s^l , this potential 
is not well-suited for MD simulations since the forces are 
discontinuous at rc [ssj . Hence, the second term is added 
in Eq. (|4]), as suggested in Ref. [1^. 

We have performed MD simulations with the BKS 
model ([T]) where we have approximated the l/r Coulomb 
interaction term by the Wolf formula, Eq. For the 
cutoff Tc in three different values have been chosen, 
namely = 6.0 A, 8.0 A, and 10. 17 A. Note that a cutoff 
of Vc ~ 10.17 A has been also used for the real space part 
of the Ewald sums in our simulations using the original 
BKS potential For a system of = 1008 particles, 
the gain in CPU time with the Wolf method is a factor of 
2 for rc = 10.17 A. Although this speed-up is not that im- 
pressive, one should keep in mind that the Wolf method 
is also well-suited for problems where Ewald sums be- 
come very inefficient, e.g. for quasi-two-dimensional ge- 
ometries, for MC simulations, or for large systems (see 
Introduction). In addition the structure of the code and 
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FIG. 2: Top: Influence of the Wolf cutoff on the Si-Si pair cor- 
relation function for large distance. Bottom: Influence of the 
Wolf cutoff on the time dependence of the self-intermediate 
scattering function for both Si (top curves) and O (bottom 
curves). Both panels are for T = 3000 K. 



the calculations of physical quantities become much sim- 
pler if the potential has only a finite range. 

Equilibrated configurations from simulations with 
Ewald sums [1^ were used as starting configurations 
for the MD simulations with the Wolf method at each 
temperature between 3000 K and 6100 K. The consid- 
ered temperatures were 3000 K, 3100 K, 3250 K, 3400 K, 
3580 K, 3760 K, 4000 K, 4300 K, 4700 K, 5200 K and 
6100 K, i.e. the same values as the ones considered in 
Ref. [3]. As in Refs. [l,!!!, the velocity form of the Ver- 
let algorithm was used to integrate the equations of 
motion with a time step of 1.6 fs. Equilibration runs were 
made by coupling the system to a stochastic heat bath, 
before performing productions runs in the microcanonical 
ensemble. During the equilibration, we detect no system- 
atic drift of the energy for rc ~ 10.17 A, which is a first 
indication that Ewald and Wolf method should yield very 
close results for this cutoff value. During microcanoni- 
cal simulations at 3000 K, the drift in the total energy is 
comparable in amplitude to the one observed when us- 
ing Ewald summations. Therefore we follow Ref. [3| and 
rescale velocities every 10^ time steps to maintain the 
total energy constant. 

We have analyzed a number of static and dynamic 
quantities (see Sec. IIIip . and we find that deviations be- 
tween Wolf and Ewald methods become more pronounced 
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when temperature decreases, and, of course, when the 
cutoff value decreases. This temperature evolution sug- 
gests that charges are more efRciently screened at high 
temperatures in amorphous silica when the system is 
more disordered. 

In Fig. [5] we present selected static and dynamic results 
at r = 3000 K for different values of the cutoff. As a dy- 
namic quantity, we show the self part of the intermediate 
scattering function (isj . 



,iq-(rj(t)- 



,(0))\ 



(5) 



where Vj (t) is the position of particle j of type a at time 
t. As a representative structural quantity we show the 
pair correlation function gapi^) (with a, = Si, O), 




r, - r, 



(6) 

where V is the total volume of the system and Na the 
number of particles of type a. In Fig. [21 the function 
gsiSi{^) is shown for r > 4 A since the Si-Si correlations 
exhibit the largest effects with respect to the cutoff Tc, 
at large distance. 

For = 6 A, deviations between Ewald and Wolf 
methods are obvious. Self-intermediate scattering func- 
tions decay faster with Wolf than with Ewald method, 
and the positions of the second and third peaks in the 
pair correlation function gsiSi('') arise at different values 
in the two methods. For Tc = 8 A, the situation is al- 
ready much better since structural relaxation occurs at 
the correct timescale, and oscillations in (7siSi(^) are in 
phase with the corresponding Ewald result. However, a 
closer inspection of the data shows that the plateau value 
in the self-intermediate scattering function is not cor- 
rectly reproduced and the amplitude of the oscillations 
in (7siSi(^) is too large at large distance. For Tc = 10.17 A, 
the agreement with the Ewald simulation is almost per- 
fect both for statics and dynamics. The height of the 
plateau in Fs{q,t) is now correct and the amplitude of 
the oscillations in gsiSi(^) are very close to those of the 
Ewald result. Looking very closely at the data, how- 
ever, a few deviations remain. The oscillations in (?siSi('') 
at large distances are still too pronounced, meaning that 
long-range order in the liquid is slightly more pronounced 
when using Wolf's method (see also the discussion con- 
cerning the results shown in Fig. [6]). The amplitude of 
the small dip in Fs{q, t) around t 1 ps is much smaller 
for = 6 A and 8 A, and remains a bit too small for 
Tc = 10.17 A, although it is only a tiny difference in the 
latter case. 

The influence of the cutoff and convergence towards 
the Ewald results is further confirmed by the behavior 
of the pressure, which will be studied in more detail in 
Sec. mil In fact, we find that the pressure represents a 
very sensitive test for the choice of Tc. For T = 3000 K, 



we find P w -3.2 GPa, -0.11 GPa, and 0.83 GPa for r,, = 
6A, 8A, and 10.17 A, respectively, while P « 0.89 GPa 
using Ewald sums for the same number of particles. 

From the results described in this section we decide 
therefore that = 10. 17 A is a good compromise be- 
tween a small cutoff which improves computational ef- 
ficiency and a very large cutoff which matches best the 
behavior observed in simulations using Ewald sums. In 
Sec. mil we shall present a more extensive set of static 
and dynamic data for Vc = 10.17 A and we will show that 
this produces a physical behavior in satisfying quantita- 
tive agreement with simulations using Ewald sums. 



C. Yukawa screening 

The replacement of the 1 /r interaction by a screened 
Coulomb interaction has been proposed in several recent 
simulation studies of various atomistic systems (see, e.g., 
Refs. [H, [13, [13, [3 and references therein). In these 
studies, screened Coulomb or Yukawa potentials have 
been considered as alternatives to the Wolf method pre- 
sented above. 

In the present work we reparametrizc the BKS poten- 
tial by replacing Vc(r) term in Eq. ([1]) by a Yukawa in- 
teraction term of the following form. 



Vyir) = Dy 



exp (—7'/ A) 



(7) 



which introduces the amplitude Dy and the screening 
length A as new parameters. Note that in Eq. ([7]) the 
same parameters Dy and A are used for Si-Si, Si-0, and 
0-0 interactions. The shape of the Yukawa interaction 
is compared to the one of the Wolf and Coulomb terms 
in Fig. m 

Our determination of Dy and A was based on previ- 
ous MD simulations for a system of 8016 particles using 
the BKS potential [J. At each temperature in the inter- 
val 6100 K< T <2750K (see Ref. ji| for the considered 
values) the parameters Dy and A are fitted such that 
the forces on each particle of the BKS configurations are 
optimally reproduced by the new potential. More pre- 
cisely, the fitting procedure was based on the following 
function. 



3N ^ a 

I a=Si,0 



(8) 



where Ff and denote the total forces on the par- 
ticle i for the original BKS potential and the new BKS 
potential modified by the screened Coulomb term ([7]), 
respectively. The (Tq represents the standard deviation 
of the BKS force distribution for Si and O particles. 
We find asi{T = 6100 K) = 4.427614 eV/A, ao{T = 
6100 K) = 3.674404 eV/A, crsi{T = 2750 K) = 3.235584 
eV/A, (7o{T = 2750 K) 2.632564 eV/A. The brackets 
(• • • ) denote an average over different samples (i.e. BKS 
configurations) . 
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FIG. 3: Dependence of the "x^ , (a), prefactor Dy, (b), and 
screening length A, (c), of the Yukawa potential (O as func- 
tions of the cutoff rc for the three indicated temperatures. 



It is important to note that if wc use the ansatz in 
Eq. ([7]) with A and Dy as the only two free parameters, 
the function could only be minimized by approaching 
the limits A = cxo and Dy = 1, thus reproducing the 
original Coulomb interaction. Therefore, we have to in- 
troduce a third parameter, namely a cutoff distance rc 
above which the potential Vy{t) is set to zero. Thus, 
we shall first fix the value of the cutoff and then de- 
termining the best set of parameters for Dy and A to 
minimize the function in Eq. ([5]). This procedure can 
then be repeated for different values of the cutoff rc, and 
for different temperatures. The outcome of this study 
will be the suggestion of an "optimal" a set of parame- 
ters (rc, Dy, A) that can be used to study the properties 



of silica at various temperatures. 

Note that the use of a cutoff in the potential ([7]) implies 
a discontinuity of the forces at rc. As already mentioned 
in the previous section, this is not useful for MD simula- 
tions since it leads to an energy drift in microcanonical 
simulation runs. Therefore, we have multiplied the forces 
by the function 



/(r) exp - 



(r - rcf 



r < rc 



(9) 



which makes the forces continuous at rc. The constant h 
was set to 2.0 A. 

The function ([5]) was minimized by means of the 
Levenberg-Marquardt algorithm [i^ which is essentially 
a conjugate gradient scheme for nonlinear fitting prob- 
lems. As shown in Fig. [3^, decreases when tempera- 
ture decreases, although the temperature dependence is 
relatively weak. It also decreases relatively quickly as a 
function of the cutoff radius rc for the Yukawa poten- 
tial. Empirically we find that the data at T = 3580 K in 
Fig. [3^ can be well described with a power law behavior, 
X^(r) oc 1/r^ '*. It is indeed expected that vanishes at 
large rc since in the limit rc — > oo, the original Coulomb 
potential should be recovered. Therefore, the screening 
length A should diverge for large values of rc while the 
amplitude Dy should approach 1, as confirmed by our 
numerical results in Figs, and c. 

Whereas the temperature dependence of Dy is rel- 
atively weak, the screening length A decreases signifi- 
cantly with decreasing temperature, which can be un- 
derstood as follows. Physically, the value of A obtained 
using the minimization procedure results from a compro- 
mise. On the one hand, a large screening length A should 
be used to recover the "bare" Coulomb potential. On the 
other hand, the introduction of a finite distance cutoff rc 
produces large errors due to incorrect charge balance [sj . 
In the present scheme, this is compensated by screening 
more strongly the Coulomb interaction by using a smaller 
screening length A. Thus the results in Fig. [5]: indicate 
that charge neutrality is best satisfied, for a fixed cutoff 
value rc, when temperature is higher, so that a larger 
screening length A can be used at high temperature. In 
turn, this suggests that charges in amorphous silica are 
screened on a smaller lengthscale when temperature is 
high, that is, when the system is more disordered. The 
same conclusion was reached in Sec. Ill Bl 

It is also remarkable that the evolution of Dy and 
A with the cutoff is not completely trivial, but exhibits 
some structure. This is best seen in the data of Fig. 
and c at r = 3580 K, for which we have taken more data 
points. This behavior can be understood from the pair 
correlation function for the SiSi correlations (Fig. 2^). 
The comparison of Fig. 0^ to Figs. and c, shows that 
the locations of shoulders or even maxima in Dy and A 
are at the locations of the minima in 5siSi('')- This makes 
sense because the minima in (7siSi(r) correspond to dis- 
tances of SiSi pairs that are relatively unlikely and so 
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FIG. 4: Influence of the cutoff Tc for the Yukawa potential on 
the Si-Si pair correlation function at T = 3250 K (a) and on 
the mean-squared displacement for the silicon atoms (b) for 
the indicated temperatures. 



nounced at the low temperature. In the diffusive regime, 
the dynamics becomes faster when decreasing the cutoff 
Tc. Thus, as in the case of the Wolf method, a change 
in long-range structural correlations is accompanied by 
a faster diffusion of the particles. In the next section, 
we present more extensive results for Tc = 10.17 A for 
the Yukawa potential. This cutoff provides a reasonable 
accuracy compared to the Ewald results, although the 
agreement cannot be made as good as the Wolf case for 
the same cutoff value. 

The simulations with the Yukawa potential (also those 
shown in Figs.|3|) were done for systems of 1152 particles 
at a total mass density of 2.37g/cm'^, i.e. the same den- 
sity that we also used for the simulations in which the 
Ewald sums and those in which the Wolf method was 
applied. Furthermore the results presented in the next 
section have been obtained considering average values for 
the Yukawa parameters, i.e. Dy = 1-07 and A = 5.649 A. 



III. DETAILED COMPARISON TO THE 
LONG-RANGE MODEL 

In this section, we present a detailed comparison 
of the static and dynamic behavior of the Wolf and 
Yukawa finite-range potentials for liquid silica and the 
one obtained with Ewald summation using parameters 
described in the previous section. 



A. Static properties 



around these distances the change in the effective screen- 
ing is weaker than for separations of SiSi pairs at which 
strong structural correlations occur. 

The pair correlation function (?sisi(r) in Fig. ex- 
hibits a similar behavior as in the case of the Wolf 
method, although the convergence towards the Ewald 
result is slower than for the Wolf method. For inter- 
particle distances r > 4A relatively strong deviations 
from the Ewald result can be observed for ~ 6.0 A. 
Note that the first peak in g<s,\<s,\(r) which is not shown 
in Fig. |3)d is also well reproduced for = 6.0 A. For 
Tc = IO.I7A the agreement with the Ewald result ex- 
tends to the second peak in gsiSi('')j while for interparti- 
cle distances r > 7A the Yukawa result is out-of-phase 
with respect to the Ewald result. 

The dynamic properties of the Yukawa potentials are 
represented in Fig. |4|3 by the mean squared displacement 
(r^(t)) for the oxygen particles (i.e. a = O) at the tem- 
peratures T = 6100 K and T = 3250 K, defined by 

(^«(i)) = ^^E|r.(i)-r.(0)py (10) 

We find similar results for Si atoms. For different val- 
ues of Tc a similar qualitative behavior at the high and 
the low temperature is seen, albeit effects are more pro- 



In Fig. Owe show the static pair correlation functions 
for Si-Si, 0-0 and Si-0 pairs at a single, low tempera- 
ture, T = 3000 K. At first sight, it is obvious that there is 
a fairly good agreement between the three sets of data. 
Looking at the data more closely, one can see that al- 
though short-range structure, < r < 5 A, is well repro- 
duced by the two finite-range potentials, small deviations 
arc present at larger distances, r > 7A, as evidenced in 
the three insets of Fig. O It is interesting that Wolf 
and Yukawa potentials produce deviations with opposite 
trends. Wolf data reveal a liquid which is more struc- 
tured at large distances than in the Ewald simulations, 
while the Yukawa potential produces less structured pair 
correlation functions. We do not have a simple physical 
explanation for this opposite tendency. We also remark 
that large distance oscillations in the Yukawa potential 
are slightly out-of-phase with the Ewald results, as is 
best seen for the Si-Si correlations (see also Fig. \^). 
But overall we conclude that the local structure of sil- 
ica is very well reproduced by the finite-range poten- 
tials. Large distance correlations are more sensitive to 
the truncation of the Coulombic interaction, but reason- 
able agreement can nevertheless be obtained for the pa- 
rameters chosen in Sec. |lll 

The opposite large distance behavior of the Wolf and 
Yukawa data are also clearly revealed when looking at 
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FIG. 5: Static pair correlation functions for Ewald, Wolf and 
Yukawa simulations for T = 3000 K (r^ = 10.17 A). The over- 
all structure of the liquid is the same for the three simulations, 
although small deviations are observed at large distances, as 
shown in the insets. 
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FIG. 6: Partial static structure factor for 0-0 and Si-Si cor- 
relations at r = 3000 K, and for the three simulations meth- 
ods. The overall structure of the fluid is correct, with small 
deviations close to the peak at 2.9 A^^. 
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the partial static structure factors, as shown in Fig. [S] 
As compared to Ewald simulations, the Wolf data cor- 
rectly reproduce S{q) for both small {q < 2.5 A~^) and 
large {q > 3.5 A~^) wavevectors, while the height of the 
peak at (7 w 2.9 A~^ is slightly too large. This is con- 
sistent with the observation of a large distance struc- 
ture seen in the pair correlation functions. Similarly, the 
Yukawa data underestimate the height of this peak, while 
its position is also slightly incorrect in the case of Si-Si 
correlations. As discussed in Sec. HIl it is important to 
recall that the amplitudes of these deviations are more 
pronounced at low temperatures when the fluid is more 
structured. 



FIG. 7: Bond angle distributions (Si-O-Si and O-Si-O) for 
two temperatures and for the three simulation methods. The 
agreement between the three models is very good at all tem- 
peratures. 

That the short-range structure of the fluid is well re- 
produced by our finite-range potentials is confirmed in 
Fig. [7] which presents Si-O-Si and 0-Si-O bond angle 
distributions for the three models. We find that at all 
temperatures the local atomic arrangements are indeed 
very similar in the three cases. We find similar results for 
Si-Si-Si angle distributions. This is consistent with the 
above observation that small differences in the pair cor- 
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FIG. 8: Temperature dependence of the pressure for the three 
potentials. 
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FIG. 10: Arrhenius plots of the self-diffusion constants for 
Ewald and Yukawa methods. 
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FIG. 9; Mean-squared displacements for the three simula- 
tions techniques at three different temperatures. At each 
temperature we show both the Si and O dynamics, the latter 
being the fastest. The agreement between Wolf and Ewald is 
excellent. The slowing down of the dynamics of the Yukawa 
potential is slightly less pronounced. 



relation functions are only seen at very large distances. 
These large distance differences do not affect the local 
tetrahedral arrangement of the atoms. 

In Fig. [S] we show the temperature dependence of the 
pressure, P. As mentioned above, we find that the pres- 
sure is very sensitive to the truncation of the Coulombic 
interaction. Nevertheless, the agreement between Wolf 
and Ewald results is very good, while somewhat larger 
deviations are observed for the Yukawa potential. Again 
the sign of the deviations is opposite for Wolf and Yukawa 
results. 



B. Dynamic properties 



In Fig. [H] we compare the mean-squared displacements 
{r^it)) for both Si and O species for several temperatures 
to the Ewald results. 

Within the statistical noise, we cannot observe any 



difference between the Ewald and Wolf data at the two 
highest temperatures. For T ~ 3000 K there is a small 
difference between the two methods on timcscales corre- 
sponding to thermal vibrations as discussed in the con- 
text of Fig. [2l A similar agreement is found for the mean- 
squared displacements in Fig. [9l although the use of log- 
arithmic scales somewhat obscures the little differences 
between the two sets of data. From Fig. [5] we also de- 
duce that the temperature evolution of the relaxation 
timescales and diffusion constants is similar using both 
Wolf and Ewald methods. It is striking that despite the 
small differences in the structure of Ewald and Wolf data, 
their dynamical behavior is in excellent agreement. This 
is not true for the Yukawa potential where larger devia- 
tions can be noticed in Fig. [9l 

The Yukawa data in Fig. [H] also imply that dynamic 
properties of the Yukawa and the Ewald models have 
a slightly different temperature evolution. This can be 
illustrated by the temperature dependence of the self- 
diffusion constants. We have calculated the self-diffusion 
constants Da via the long-time limit of the mean squared 
displacements using the Einstein relation: 



rl{t)) 



Da = lim 



(11) 



Figure [TO] displays the diffusion constants in an Arrhe- 
nius plot. Note that we do not plot the Wolf results 
in this figure, because there are only minor differences 
between Wolf and Ewald data in the full temperature 
range. Whereas at the highest temperature, T ~ 6100 K, 
Yukawa and Ewald methods yield the same values for 
Da within the statistical accuracy, at lower temperatures 
the Yukawa values are a factor of 2 to 3 higher than the 
"exact" Ewald values. Note that at low temperatures 
the Yukawa result for D^i coincides with the Ewald re- 
sult for Do- This is certainly a coincidence but it il- 
lustrates that qualitatively the temperature dependence 
as obtained from the Yukawa potential is very similar 
to that of the "exact" calculation with Ewald sums. In 
particular, from both methods similar activation ener- 
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FIG. 11: Vibrational density of states for the three potentials. 

gies are expected for the temperature dependence of Dq, 
in the low-temperature regime. 

Finally we show in Fig. [Tl] that the vibrational den- 
sity of states 5(1/) is well reproduced by the Wolf method 
and, on a qualitative level, also by the Yukawa method. 
In order to compute g{i') we quenched 10 independent 
liquid samples with an infinite cooling rate, i.e. a steep- 
est descent, down to T = BOOK followed by an annealing 
during 80 ps. We then measured the velocity autocorre- 
lation function and obtained g{iy) from its Fourier trans- 
form [4^. The liquid starting configurations were well- 
equilibrated samples at T = 2715 K for the Ewald as well 
as the Wolf case and at T = 3000 K for the Yukawa case. 
As we see in Fig. [11] Ewald and Wolf results for 5(1/) 
are in very good agreement, while the Yukawa method 
yields to slight deviations from the Ewald calculations, in 
particular at high frequencies {ly > 30 THz) where gi^v) 
exhibits a double-peak structure. The latter peaks are 
due to inter-tetrahedral stretching modes. In the past it 
has been shown that cooling rate effects affect especially 
the amplitude of the high-frequency peaks 0,0. Thus, 
part of the deviations seen in the Yukawa results might 
be due to the different cooling history that we used for 
the Yukawa case. Wc also mention that in the frequency 
range 5 THz< 1/ <25 THz the vibrational dynamics in the 
BKS model is quite unrealistic, as has been revealed by 



the comparison to Car-Parrinello MD simulations [l3l |. 
The results of the present study suggest that an improve- 
ment of the BKS model with respect to vibrational prop- 
erties can be also achieved by a finite-range potential. 



IV. CONCLUSION 

Extensive MD simulations have been used to study 
finite-range approximations of the BKS model for sil- 
ica. Wc have demonstrated that both the Wolf method 
and a screened Coulomb (Yukawa) potential can be used 
to approximate the Coulomb interactions in BKS silica. 
In order to obtain quantitative agreement with Ewald 
sums, a cutoff of about 10 A has to be used for the 
Wolf potential. Such a potential cannot be classified 
as a short-range potential since its spatial extent corre- 
sponds to about three connected Si04 tetrahedra. How- 
ever, due to its finite range, the Wolf potential can be 
efficiently applied to problems where Ewald sums be- 
come inefficient, such as quasi-two-dimensional geome- 
tries, MC simulations 0, or MD simulations with a 
very large number of particles. In the case of the Yukawa 
potential larger cutoffs have to be used to yield quanti- 
tative agreement with Ewald results. However, we have 
to keep in mind that we have simply reparametrized the 
Coulomb interaction term of the BKS model. An inter- 
esting possibility would be to perform a more complete 
reparametrization of the BKS potential, so that a smaller 
cutoff value could be used, making simulations even more 
efficient. 
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